clear;

%Calculate change of a intensity profile after transimitting a multiple thin film system 
%Written by Feiran Huang on 05/31/99. Update JW Haus 11/16/99 
% Updated and corrected by JW Haus on 8/21/2002 
%------------------------- 
%given parameters for film system, n,k,d of various wavelength 
%given expression of input intensity profile in time domain and center wavelength  
v=[500 2100 0 1];  
%film information 
d1=1000/4/2; 
d2=1000/4/1.;
% number of periods
m=10;
%  angle of incidence 
fi0=0/180*pi; 
n00=[1. 1. 1. 1.  1.];
lambda=v(1):.5:v(2);  
ld0=[180 500 600 800 6100]; 
n10=[2 2 2 2 2];
k10=[0.0 0.0 0.0 0.0 0.0];
%k10=[0.01 0.01 0.01 0.01 0.01]; 
n20=[1.0 1.00 1.00 1.00 1.00]; 
k20=[0.0 0.0 0.0 0.0 0.0];
%k10=[-0.01 -0.01 -0.01,-0.01 -0.01]; 
n30=[1. 1. 1. 1. 1.]; 
%n30=[1. 1. 1. 1. 1.];

n0=interp1(ld0,n00,lambda);
n1=interp1(ld0,n10,lambda); 
k1=interp1(ld0,k10,lambda); 
n2=interp1(ld0,n20,lambda); 
k2=interp1(ld0,k20,lambda); 
n3=interp1(ld0,n30,lambda);  
[Tp,Rp,Ts,Rs]=f_2mlyr(n0,n1,k1,n2,k2,n3,d1,d2,m,fi0,lambda);  
figure(1);
plot(lambda, Tp,'b',lambda, Ts,'r');
ylabel('Transmittance'); 
xlabel('Wavelength (nm)'); grid on; axis(v);
figure(2);
plot(lambda, Rp,'b',lambda, Rs,'r');
%plot(lambda, Tp,lambda,Rp,lambda, Rp+Tp); 
%data0(:,1)=lambda.'; 
%data0(:,2)=pTransmittance.'; 
%save data01 data0 -ascii; 
ylabel('Reflectance'); 
xlabel('Wavelength (nm)'); grid on; axis(v);

% field amplitude caculation
lambdae=1575;%1312;%1413;%1575;
n0e=interp1(ld0,n00,lambdae);
n1e=interp1(ld0,n10,lambdae); 
k1e=interp1(ld0,k10,lambdae); 
n2e=interp1(ld0,n20,lambdae); 
k2e=interp1(ld0,k20,lambdae); 
n3e=interp1(ld0,n30,lambdae);  
[Ees,Eep,fi1,fi2,fi3]=f_2mlyrAmp(n0e,n1e,k1e,n2e,k2e,n3e,d1,d2,m,fi0,lambdae);
count=1;
dz1=d1/10;
dz2=d2/10;
z(1)=0;  % input boundary
Eys(1)=Ees(1,2)+Ees(2,2);
Exp(1)=Eep(1,2)+Eep(2,2);
count=count+1;

for jj=1:m          %check for the intial value 1 or 2
   
   % layer 1
   
  for kk=0:9
Eys(count)=(((Ees(1,2*jj)*exp(-i*2*pi*(n1e-i*k1e)*cos(fi1)/lambdae*kk*dz1))+(Ees(2,2*jj)*exp(i*2*pi*(n1e-i*k1e)*cos(fi1)/lambdae*kk*dz1))));
Exp(count)=(((Eep(1,2*jj)*exp(-i*2*pi*(n1e-i*k1e)*cos(fi1)/lambdae*kk*dz1))+(Eep(2,2*jj)*exp(i*2*pi*(n1e-i*k1e)*cos(fi1)/lambdae*kk*dz1)))*cos(fi1));
 z(count)=z(count-1)+dz1;
 count=(count+1);
end

% layer 2
   

for kk=0:9
Eys(count)=(((Ees(1,2*jj+1)*exp(-i*2*pi*(n2e-i*k2e)*cos(fi2)/lambdae*kk*dz2))+(Ees(2,2*jj+1)*exp(i*2*pi*(n2e-i*k2e)*cos(fi2)/lambdae*kk*dz2))));
Exp(count)=(((Eep(1,2*jj+1)*exp(-i*2*pi*(n2e-i*k2e)*cos(fi2)/lambdae*kk*dz2))+(Eep(2,2*jj+1)*exp(i*2*pi*(n2e-i*k2e)*cos(fi2)/lambdae*kk*dz2))).*cos(fi2));
 z(count)=z(count-1)+dz2;
 count=(count+1);
  end
end
   
   
   figure(3)
   plot(z,abs(Exp).^2)
ylabel('|E|^2 P-polarization'); 
xlabel('Z (nm)');
figure(4)
   plot(z,abs(Eys).^2)
ylabel('|E|^2 S-polarization'); 
xlabel('Z (nm)');
